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We study the multicritical behavior for the semimetal-insulator transitions on graphene’s honey¬ 
comb lattice using the Gross-Neveu-Yukawa effective theory with two order parameters: the SO(3) 
(Heisenberg) order parameter describes the antiferromagnetic transition, and the Z 2 (Ising) order 
parameter describes the transition to a staggered density state. Their coupling induces multicritical 
behavior which determines the structure of the phase diagram close to the multicritical point. De¬ 
pending on the number of fermion flavors Nj and working in the perturbative regime in vicinity of 
three (spatial) dimensions, we observe first order or continuous phase transitions at the multicritical 
point. For the graphene case of Nj = 2 and within our low order approximation, the phase diagram 
displays a tetracritical structure. 


I. INTRODUCTION 


Graphene [T] features a number of unconventional elec¬ 
tronic properties which to a large extent can be captured 
by means of a simple tight-binding description of the elec¬ 
trons on the honeycomb lattice |2]. The experimental 
findings such as, for example, the half-integer quantum 
Hall effect [5] or the Klein paradox [3] suggest that this 
single-particle picture provides a very good starting point 
for the theoretical description of the material. Electron- 
electron interactions therefore are often expected to play 
only a quantitative and not a qualitative role. In fact, on 
the honeycomb lattice at charge neutrality the density 
of states vanishes linearly for energies close to the Fermi 
level, and an interaction-induced transition toward an 
ordered - possibly Mott insulating - state appears only 
when a minimal critical value of the interaction strength 
is exceeded [SHE]. Depending on the interaction profile, 
i.e. the precise ratios of onsite, nearest-neighbor, and fur¬ 
ther interaction parameters, a great variety of different 
spontaneously broken symmetries can occur [me]. The 
symmetry breaking patterns, most prominently, include 
chiral symmetry breaking phases, such as the antiferro¬ 
magnetic spin density wave (SDW) or a charge density 
wave (CDW) [TBHIS]. However, more exotic states of 
matter such as quantum Hall and the quantum spin Hall 
phase [SHII], or the existence of a quantum spin liquid 
have been discussed [501111]. 

Recent ab initio calculations [HHIO], on the other hand, 
find values for the interaction parameters of graphene 
which seem to be close to a transition toward a sponta¬ 
neously symmetry-broken phase. This motivates a con¬ 
sideration of circumstances under which graphene could 
be driven through a quantum phase transition by tuning 
of external parameters. One possibility could be a uni¬ 
form mechanical stretching of the graphene sheet to in¬ 
crease the ratio of onsite interaction to nearest-neighbor 
hopping parameter. The ratio between the onsite and 
the nearest-neighbor interaction might also be in prin¬ 
ciple tuned by proximity of a substrate that induces a 
screening of the Coulomb interaction. 

Here, we consider a situation with both a large onsite 


interaction, and a sizable nearest-neighbor interaction. 
Separately, these interactions would trigger phase transi¬ 
tions toward the ordered SDW and the CDW states, re¬ 
spectively. The experimental findings together with the 
ab initio parameters suggest that graphene is close, but 
somewhat below the critical values for the formation of 
one of these ordered states, cf. Fig. [T] 

The separate critical behavior of these two orders for 
Dirac electrons on graphene’s honeycomb lattice has been 
investigated in many studies; see Refs. |T6l |T9l |24l [25l for 
example. On the other hand, with both interaction pa¬ 
rameters being close to criticality, the system is close to 
the point where the transition lines into the SDW and 
the CDW phases meet, and multicritical behavior results. 
The problem of competing order parameters is a com¬ 
plex many-body problem, and multicritical behavior is 
important to many different fields of physics. A promi¬ 
nent example arises in the field of correlated electrons, 
e.g. in the study of high-temperature superconductors 
and related complex materials. In graphene, multicriti¬ 
cal behavior for the case of broken spin rotation symme- 



FIG. 1: (a) Schematic phase diagram of the extended Hub- 
hard model on the honeycomb lattice with onsite interaction 
U and nearest-neighbor interaction V. Neutral suspended 
graphene is found to be in the semi-metallic state indicated 
by the star. The neighborhood of the multicritical point (grey 
shaded area) may be governed by a (b) tetra-critical or (c) 
bicritical structure or by a (d) first-order multicritical point. 
Solid lines denote second order and dashed lines first order 
phase transitions. 







2 


try has been studied [ 55 ]. At the same time the study 
of low-dimensional Dirac electrons and their competing 
interactions can shed light on related structures in the 
phase diagram of quantum chromodynamics. 

In this paper, we study the multicritical behavior of an 
effective theory for the semimetal-insulator transitions on 
graphene’s honeycomb lattice, which assumes the form of 
a Gross-Neveu-Yukawa field theory with two coupled or¬ 
der parameters. We consider a generalized theory of this 
type with an arbitrary number of fermion flavors Nf. 
The SO( 3 ) order parameter describes the antiferromag¬ 
netic transition and the Z2 order parameter describes 
the transition to a staggered density state. Their cou¬ 
pling induces multicritical behavior which determines the 
structure of the phase diagram close to the multicritical 
point | 27 j . A first order transition at the critical point 
is possible, or if the transition is continuous, potential 
bi- or tetracritical behavior decides over the existence of 
a mixed phase in the symmetry-broken regime (see in¬ 
set in Fig. [^. We set up an e expansion close to three 
(spatial) plus one (time) dimensions to study which one 
of these three possibilities is realized for electrons on the 
honeycomb lattice, and, in particular, investigate the de¬ 
pendence of the result on the number of fermion flavors 
Nf. For small flavor number, to the first order in e, we 
find that the multicritical behavior is governed by the 
decoupled fixed point that consists of the chiral Heisen¬ 
berg and the Ising fixed points from the two separate 
descriptions of the Mott insulator transitions. This im¬ 
plies tetracritical behavior and a decoupled mixed phase 
of the two order parameters. For intermediate fermion 
flavor numbers, all transitions are of first order. Finally, 
for a large number of flavors the phase diagram displays 
a bicritical structure with a first order transition between 
the CDW and SDW state. 

The remainder of the paper is organized as follows. In 
the next section we present the effective model with two 
dynamical order parameters. We then compute the fixed- 
point structure as a function oi Nf within first-order e 
expansion in Sec. |TTT| In Sec.jTVjwe discuss the resulting 
phase diagram and give concluding remarks. 

II. EXTENDED HUBBARD MODEL ON THE 
HONEYCOMB LATTICE 

To describe the behavior of electrons on the honey¬ 
comb lattice we start with the tight-binding Hamiltonian 
supplemented by the interaction terms, H — Hq + H-mt, 
with 

Ho =-t ul{R)vs{R + 6i)+ h.c., (1) 

R.Si.s 

Hint — U ^ F F ^ t (5) 

where Ug and Vg are the electron annihilation operators 
at the two triangular sublattices of the honeycomb lat¬ 


tice with spin projection s = f, j, and R denotes the 
sites of one triangular sublattice. The 6 i are the vec¬ 
tors to the three nearest neighbors on the second sub¬ 
lattice. Explicitly, the position vectors of the bipar¬ 
tite lattice are spanned by Rf = a (■\/ 3 / 2 , — 1 / 2 ) and 
R2 = a (0,1), where a is the lattice spacing which in 
the following is set to a = 1 . The second sublattice is 
generated by i? -|- with the three nearest-neighbor vec¬ 
tors = ( 1 /( 273 ), 1/2), ^ = (1/(273),-1/2), and 
( 5 |’ = (1/73,0). Then, for the non-interacting part of 
the Hamiltonian, the energy spectrum is described by 
two bands eg = ±t| exp(fc • ( 5 i)| which are linear and 
isotropic close to zero energy near the K, K' points at the 
border of the Brillouin zone, where = ( 27 r/ 73 , 27 r/ 3 ) 
and K' = —K. Employing only Fourier-modes near the 
AT, K' points the continuum low-energy effective action 
corresponding to Ho can be written as [ 5 ] 

/.l/T 

Sf= dTdx^~^ [^' (I2 0 7 At) ( 3 ) 

Jo 

where D is the space-time dimension and the spin-^ 
electrons are described by the 8-component Dirac field 
xj/ = (xh^,d>^)^ with its conjugate ^ G) 7o) in 

2 < D < 4 . Further, we have 5 ^ = ( 9 i-,V) and the 
Clifford algebra {7^,71/} = 25 ^i, with (4 x 4)-7 matrices 
and we assume the summation convention over repeated 
indices, and z/ = 0 ,..., D— 1 . Explicitly, in ( 2 -|-l)D we 
use a representation where 70 = I2 <8 7i = 0 Oj,, 

72 = 12 002;. With these definitions the relation between 
the Grassmann fields u, v and can be given as [ 5 | 

Tt(F,r) = tJ2J 

vl{K + q,uJn),ul{-K + q,uJn),vl{-K + q,ujn)\. ( 4 ) 

The reference frame is chosen to be such that qx = 
q ■ K/\K\ and qy = {K x q) x K/\K\‘^ and further we 
set h = ks = vp = 1 . We also dehne the two additional 
4 x 4 matrices that anticommute with all 7^, namely 
^0 = cTx® ay and 75 = ay ay. Then 735 = -i'folo com¬ 
mutes with all 7^ and anticommutes with 73 and 75. In 
the following, we will also consider the generalization to 
arbitrary number of Dirac points in the spectrum. For¬ 
mally, this can be done by replacing Us{zLK + g,a;„) 1—>■ 
Ug{±Kt + g, Wn) and Vg{±K -|- g,a;„) >-)■ Vg{±K^ + q,uin) 
with Ki the position of the Dirac points, i = 1 ,..., Nf. 
We will refer to Nf as the number of fermion flavors with 
Nf = 2 being the physical case. For general Nf, the gen¬ 
eralization of Eq. (^, thus renders the fermion field to 
have 2Nf spinorial components for each spin direction. 

Order Parameters and Interactions There are several 
different order parameters that induce various symmetry 
breaking patterns, e.g., time-reversal symmetry (TRS) 
breaking and chiral symmetry (CS) breaking. In this 
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work, we will consider CS breaking while TRS is pre¬ 
served. The CS breaking order parameter can be written 
as 


$= (X,<^ = l2iv,)«')) . (5) 

The Ising field x is a spin singlet, and corresponds to a 
staggered density, or a charge density wave state, x 
ulus — vlvg, which can be triggered by a large nearest- 
neighbor density-density interaction. The Heisenberg 
field (/) is a triplet and corresponds to a staggered magne¬ 
tization or an antiferromagnetic spin density wave state, 
(p ~ ulass'Us + vlass'Vs, and it is triggered by a strong 
onsite interaction. The order parameters which appear 
in the form of fermion bilinears can be promoted to be 
dynamical fields corresponding to a bosonic action 


Sb= / drdx 


•D-i [1,,/ q 2 , ,1 r / c,2 


■ ( 6 ) 


Finally, there are also Yukawa terms that couple bosonic 
and fermionic degrees of freedom 


Sy= Jdrdx^ ^ ® 

-I- g^$- §((T 0 l2Af/)'I' 


(7) 


The complete action then is given hy S = Sp + Sb+ Sy- 
Its form is dictated by the spin-rotational, time-reversal 
and sublattice-exchange symmetries. This effective the¬ 
ory is explicitly relativistic, cf. Refs. [HI [19]. The 
Coulomb repulsion does not appear explicitly, while it 
may tune the transition entering through the nearest- 
neighbor repulsion. Its long-range tail has been shown to 
be an irrelevant perturbation at the critical points [161 

I25j. 


III. e EXPANSION 

For the above action S and at zero temperature we 
can calculate the renormalization group equations in the 
Wilsonian scheme by simultaneously integrating out the 
fermionic as well as the bosonic modes within the narrow 
momentum shell A/& < + (f) < A. At one-loop order 

in Z) = 4 — e dimensions we find the following equations 
for the two Yukawa couplings 

^ = e4-(3 + 2iV/)4-95^4, ( 8 ) 

dq"^ 

^ - (1 + 2Nf)gl - iglgl , (9) 

where we have rescaled gf —>■ gf /A’^). Addition¬ 
ally, in the e expansion we obtain RG flow equations 
for the bosonic masses and the bosonic couplings 


Ax, A,^, Ax 0 that will be displayed below. The fixed points 
of this set of equations and their surrounding will deter¬ 
mine the properties of the multicritical point. 

The flow in the g‘^-g\ plane decouples from the bosonic 
flow equations, and analytical solutions of the zeros of the 
Yukawa-coupling beta functions {, 5 ^ 2 , / 3 g 2 } can be readily 
displayed. We find the values 


A-.g 
B: <7 

C-.g 


D: 5 


= 0 , gY = 0 , 


2,* ^ Q 2.* ^ _ _ 

X ^ y<p 1 -h 2 A/ ’ 


2,* p. 

iX = 0 , 


2 ,* 

X 


3 + 2Nf' 


2 ,* _ 2 
X “ 


{Nf - 4)e 


Nf + 2Nf-6' 


9^ 


( 10 ) 

( 11 ) 


( 12 ) 


\Nf 


Nj + 2Nf-Q' 


(13) 


We illustrate the flow and fixed-point structure in the 
50-l?x different values of Nf in Fig. Let us 

emphasize, however, that the zeros A, B, C, and D, in 
order to represent true fixed points of the full system, 
need to be supplemented with suitable (and physically 
admissible) fixed-point values in the bosonic sector. In 
terms of the rescaled bosonic couplings {A^^, A^, A^^} —)■ 
{Ax/( 87 r^A'^), A,^/( 87 r^A'^), Ax 0 /( 87 r^A'^)} the /3 functions 
in the bosonic sector read 


din 6 
dX(p 
din 6 
dA(^x 
din 6 


= eAx - 36A^ - - iNfX^g^. + Afyg'^ , 

= eX^ — 44 A 0 — A^^ — ANfXfj^g^ + Nyg'^ , 
= eA^x ~ SA^^ — 2 OA 0 A 0 X ~ 12 AxA,^x 


f^<t>x94> ‘^^f^<t>x9x 9x9(1) ■ 


(14) 

(15) 

(16) 


For a fixed point A, B, C or D to be physically ad¬ 
missible it has to satisfy a number of conditions. First, 
the square of the Yukawa couplings has to be positive, 
gf > 0,i G {<(), x}- Also we need X^,X^ > 0 to have 
an action that is bounded from below. This is accompa¬ 
nied by the condition that A^x is not too negative, i.e. if 
A^x < 0 we need dAxA^ — A^^ > 0. Further, we will have 
at least two relevant directions related to the mass pa¬ 
rameters of the bosons x and (p. They correspond to the 
two tuning parameters in the phase diagram in Fig. 
This implies that the physical fixed point may not have 
more than these two relevant directions, in order to be 
accessible without tuning a third microscopic parameter. 
Such fixed points will be called here “stable”. We hence 
classify the hxed points according to their number of rel¬ 
evant directions, i.e., their number of positive eigenvalues 
of the stability matrix. For the fully Gaussian fixed point 
A, we hnd that the Yukawa couplings already provide two 
additional relevant directions with eigenvalues e for any 
Nf. Thus in the following it will be discarded in our 
discussion. 

Fixed point B leads to the least number of relevant 
directions in the full system when supplemented with the 
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Nf < 4 Nf = 4 






FIG. 2: RG flow and fixed-point structure in the Yukawa-coupling subsector for different values of Nf below, at, and above the 
critical value of Nf = 4. For Nf < 4 fixed point B is stable within this subsector, while fixed point D is located in the unphysical 
domain < 0 (top left). At Nf = 4, D becomes physical and collides with B (top right), and subsequently exchanges stability 
with the latter for Nf > 4 (bottom left). For Nf — >■ oo the stable fixed point D is located at = g^ (bottom right). Note 
that A, B, G, and D may be true (and possibly stable) fixed points of the full system only if suitable corresponding zeros of 
the bosonic beta functions can be found, see text. 


bosonic couplings 


a; 


36’ 


a;^ = o, 


1 - 2iV/ -b y/l + 4Nf{43 + Nf) 
88{2Nf + 1) 


(17) 

(18) 


which, in fact, means that the order parameters (j) and x 
decouple at B, with the Heisenberg field at its fermionic 
chiral fixed point, and the Ising field at the purely 
bosonic (standard) Ising fixed point. We thus denote 
fixed point B also as “chiral Heisenberg plus Ising” 
(cH-l-I) fixed point. Its critical exponents are 


Oi 

O 2 

O 3 


= 2 - 

3 ’ 

5 -b 347V/ -b Sv'l + 47V/(43 -b Nf) 

— 2 — ^ 7^ c, 


22i2Nf -b 1 ) 


Nf-A 


(19) 

( 20 ) 
( 21 ) 


For the physical situation with Nf = 2 we obtain 
{ 61 , 02 , 63 } = {2 - |,2 - 1^,-^}, rendering the cH-bl 
fixed point stable as 63 < 0. The situation changes upon 
increasing the value oi Nf, cf. Fig. For Nf = 4, B col¬ 
lides with another fixed point D, with which it exchanges 
stability for Nf > 4. For Nf > A the cH-bl fixed point 
thus becomes unstable. The scenario is similar to the 
well-known situation in the purely bosonic 0{N) mod¬ 
els, in which the Wilson-Fisher fixed point approaches 
the Gaussian fixed point when the dimension is increased 
towards the upper critical dimension, and subsequently 
exchanges stability with the latter. 

Fixed point C, when completed with the bosonic fixed- 
point values 


Ax = 


3-2Nf + y^9 + ANf{88 + Nf) 
72{8 + 2Nf) ' 


( 22 ) 


X* — 


'i> 44 ’ 


a;^ = 0 . 


(23) 
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TABLE I: The largest three critical exponents and anomalous 
dimensions for fixed points B and C in first order e expansion 
for the physical choice Nf = 2. Here, fixed point B (“chiral 
Heisenberg plus Ising” fixed point) is stable and thus governs 
the multicritical behavior. In the decoupled bosonic sector 
the anomalous dimension vanishes as it is of 0{e^). 



01 02 

03 

Vf, 

Vx 

ili> 

B 

(cH-l-I) 

2 - 2 - ||e 

-1^ 


0 


C 

(cI+H) 

2-^f 2-^f 


0 




could then, within our nomenclature, be termed “chiral 
Ising plus Heisenberg” (cI+H) fixed point. It has three 
relevant directions for all admissible choices ot Nf and e. 


01 

02 


= 2 -^, 

11 ’ 

3 + lONf + y^9 + 4fV/(33 + iV/) 
“ ^ 6(3 + 2Nf) 


03 


2Nf 

3 + 2iV/ ^ ’ 


(24) 

(25) 

(26) 


and is thus never stable. Its critical exponents and 
anomalous dimensions for the physical case Nf = 2 are 
listed together with the those of fixed point B in Table |l] 


For the mixed fixed point D the computation of the 
bosonic fixed-point values A* is slightly more involved, 
but can readily be done numerically. The bosonic fixed- 
point couplings lead to a quite intriguing behavior as a 
function of Nf, as displayed in Fig. First, we observe 
that there are two ranges of Nf where no real fixed-point 
values in the bosonic sector can be found for the given 
values of and 5 *^ [Eq. @]. The first range covers 
small Nf < 3.8 and the second range is Nf G [4.7,15.7]. 
These intervals are indicated in Fig.j^by the gray-shaded 
areas. In the narrow range Nf G [3.8,4.7] the solution for 
fixed point D is physically admissible only for Nf > 4, 
since < 0 for iVy < 4. For Nf = 4 the fixed points 
B and D collide and exchange stability, so that for 4.7 > 
Nf > 4 the fixed point D is stable. For 15.7 < Nf < 16.6, 
while corresponding real fixed-point values can now be 
found, the solution for fixed point D remains unphysical, 
since A;^ < 0 (see Fig. [3). Finally, for Nf > 16.6 the fixed 
point D constitutes a ^ysically admissible solution with 
negative exponent 03 , i.e., it is the stable fixed point. 
We have checked that between Nf = 4 and iVy = 16 no 
other real and stable solution exists in our approximation 
scheme. In the limit of large Nf, we find 


2,*_ 2,* _ £ 

-9<j> - ^ ’ 


(27) 


and we obtain a stable fixed point in agreement with 
the large-calculation [8]. It has the bosonic coupling 
coordinates 


A*=A 


* 


e _ 3e 


(28) 



FIG. 3: Upper panel: Coordinates of the bosonic couplings at 
the fixed point D as a function of Nf. The gray areas indicate 
the ranges of Nf where fixed point D does not exist and no 
other stable fixed point exists. Lower panel: Third critical 
exponent as function of Nf for the two possible stable fixed 
points B and D. 


IV. DISCUSSION 

An important quantity to classify the critical behavior 
of a multicritical point is the parameter A, 

A = 4A;,A^ - A|^, (29) 

whose sign determines the nature of the multicritical 
point, i.e., whether it is bicritical or tetracritical. When 
A < 0 the transition line between the two ordered phases 
is expected to be of first order, whereas for A > 0 coex¬ 
istence of the two orders is expected with second order 
transitions between the four different regimes. This con¬ 
clusion is based on a consideration of the saddle-point 
approximation of the free energy |291 130) . and we ex¬ 
pect that the presence of Dirac fermions does not affect 
it, since the argument relies only on the consideration 
of the boson effective potential. We identify three dif¬ 
ferent regimes if we classify the nature of the multicrit¬ 
ical point in terms of A, cf. Fig. For small fermion 
flavor number Nf < 4.64, including the graphene case 
Nf = 2, we find a positive A at the stable fixed point. 
Here, tetracritical behavior dominates the phase transi¬ 
tions and a mixed phase appears with a coexistence of 
SDW and CDW order. In a large range of this regime 
1 < Nf < 4, the cH-|-I fixed point B with vanishing is 
stable indicating that the properties at the multicritical 
point are dominated by the chiral Heisenberg universal¬ 
ity class. This universality class also describes the tran- 
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FIG. 4: Value of the parameter A as a function of Nf for the 
stable fixed point, which governs the multicritical behavior 
in the phase diagram of Fig. j^a). For Nf < 4 the mul¬ 
ticritical point is given by the chiral Heisenberg plus Ising 
fixed point (B) and it has positive A > 0, corresponding to 
tetracritical behavior with a mixed CDW and SDW phase, 
as indicated in Fig. Bb). For 4: < Nf < 4.64 the multi¬ 
critical point is still tetracritical, but its long-range behavior 
defines a new universality class governed by fixed point D. 
For 4.64 < Nf < 4.7 the multicritical point is bicritical (see 
Fig- Be)). For 4.7 < Nf < 16.6 there is no physically ad¬ 
missible stable fixed point in the system and we expect the 
multicritical point to become discontinuous. For Nf > 16.6 
the multicritical behavior is again bicritical. 


sition from the semimetal to the antiferromagnetic state 
for increasing nearest-neighbor interaction. The stability 
suggests the continuation of this behavior up to the mul¬ 
ticritical point. On the bosonic side of the transition the 
order parameters decouple, whereas they interact with 
each other if fixed point D becomes stable. The sta¬ 
bility of fixed point D is accompanied by a negative A 
for Nf G (4.64,4.7) and Nf > 16.6. In this regime the 
phase diagram exhibits bicritical behavior with continu¬ 
ous transitions between SM and SDW or CDW, respec¬ 
tively, and a first order line between SDW and CDW. 
The large-fV/ behavior is thus consistent with the pre¬ 
vious i/Nf approaches [5]. The intermediate regime 
Nf G (4.7,16.6) is characterized by the absence of any 
stable fixed point. Here, the second order lines of the 
separate transitions from SM to SDW and CDW end in 
a discontinuous point and the transition from SDW to 


CDW is also of first order. 

In conclusion, we presented an effective description 
of the multicritical point between the semimetallic, the 
CDW and the SDW state on the honeycomb lattice that 
becomes exact near three spatial dimensions. A rather 
complex picture emerges as a function of the fermion fla¬ 
vor number, in which the graphene case appears domi¬ 
nated by the chiral Heisenberg universality class, at least 
to the present order of the e expansion. A tetracriti¬ 
cal phase diagram with a coexistence phase of SDW and 
CDW is the consequence. In purely bosonic multicritical 
systems [23 130] , as well as in the separate Gross-Neveu- 
Yukawa models |19j . the leading-order e expansion cap¬ 
tures well the qualitative behavior of the systems. We 
expect this to hold for our model as well and, at the 
very least, our results should qualitatively correct de¬ 
scribe the Nf dependence of the fixed-point structure in 
two spatial dimensions as appropriate for the description 
of graphene. However, the bare numbers of the critical 
Nf values at which the multicritical behavior of the sys¬ 
tem changes can possibly be reduced substantially if one 
went beyond the leading-order e expansion 133133 and 
more elaborate investigations may be needed to settle 
the true nature of the multicritical point in graphene’s 
phase diagram. A promising complimentary approach 
to the critical Nf values could be provided by a func¬ 
tional renormalization group study |19[ I33j which allows 
to work directly in two spatial dimensions. This is under 
way [54] . 
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